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v/-^ I The mass function of 51 Massive Compact Objects (MACHOs) detected toward 

^ ■ the Galactic bulge is statistically estimated from Einstein ring crossing times t^. For 

^ ■ a Gaussian mass function, the best fitting parameters are (log(m/M0)) = —1.12 and 

criog(m/M0) = 0.57. If the mass function follows a power-law distribution, the best 
2 ! fitting mass cut-off and power are rricut = 0.04 Mq and p = —2.1. Both mass spectra 

^ \ indicate that a significant fraction of events are caused by MACHOs in substellar mass 

■ range. Both best determined mass functions are compared with that obtained from 

■ Hubble Space Telescope (HST) observations. The power law is marginally favored 
^ . (2.3o") over the Gaussian mass function, and significantly over HST mass function 

(5.5(t), indicating that MACHOs have a different mass function from stars in the 
solar neighborhood. In addition, the fact that all the models have very poor fits to 
the longest four events with tf. > 70 days remains a puzzle. 
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1 Introduction 



Major efforts to detect MACHOs by observing microlensing events of source stars 
located in tlie Galactic bulge have been carried out by the OGLE (Udalski et al. 1994) 
and MACHO (Alcock et al. 1995a) groups. Current and prospective observations can 
constrain the disk/bulge normalization (Stanek 1995) and the mass distribution of 
the bulge (Kiraga & Paczyhski 1994; Evans 1994; Han & Gould 1995b). 

However, it is very difficult to obtain information about the physical parameters 
of the individual lenses. This is because the only measurable quantity from current 
observations, the Einstein ring crossing time t^, depends on a combination of physical 
parameters of the individual lenses. The Einstein ring crossing time is related to the 
physical parameters of the lenses by 
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where m and v are the mass and the transverse speed of the lens, is the Einstein 
ring radius, and Dgi, D^s, and Dig are the distances between the observer, lens, and 
source. 

Several ideas have been proposed to break the degeneracy of t^. The distance 
to a lens can be measured when the lensing star is bright enough to be detected 
(Kamionkowski 1994; Buchalter, Kamionkowski, & Rich 1995). The ambitious and 
promising idea of measuring MACHO parallaxes from a satellite would clarify the 
dynamical motions of MACHOs (Gould 1994b, 1995; Han & Gould 1994a). The MA- 
CHO proper motion, fi = v/Dqi, can be measured photometrically (Gould 1994a; Ne- 
miroff & Wickramasinghe 1994; Witt & Mao 1994; Witt 1995; Loeb & Sasselov 1995; 
Gould & Welch 1996) and spectroscopically (Maoz & Gould 1994). Han, Narayanan, 
& Gould (1996) also proposed a method to measure /i by using lunar occultation. 
When the information from proper motion and parallax measurements are combined, 
the degeneracy of can be completely broken, yielding m, and Dd- However, 
measuring the distance or proper motion of the lens star is possible only under very 
restricted conditions. Most lenses are hkely to be too faint to be detected. A proper 
motion is measurable only when a lens crosses over or very close to the face of the 
source star. The Lunar occultation method is applicable only to a small number of 
events becaiisc the hmar path is restricted. Although satellite-based parallaxes are 
promising, actual measurements are many years off. 

However, one can still obtain much information about the mass spectrum of MA- 
CHOs from the data on time scales that are provided by current observations: MA- 
CHO (Alcock, et al. 1995b) and OGLE (Udalski et al. 1994) have detected > 50 



events during 1993 bulge season. De Riijula, Jetzcr, & Masso (1991) have developed 
a method of "mass moments" to analyze MACHO masses and Jetzer (1994) has ap- 
plied this to measure the mean mass (0.28 Mq) of bulge MACHOs using 4 OGLE 
events. 

In this paper, we use maximum likelihood to estimate the mass function of mi- 
crolenses which are detected toward the Galactic bulge by combining the latest te 
data from MACHO and OGLE with plausible models for the velocity and spatial 
distributions of the lenses and sources. Wc test both Gaussian and power-law mass 
functions, and find (log(m/M0)) = —1.12 and (Tiog{m/MQ) = 0.57 for a Gaussian mass 
function, and mcut = 0.04 Mq and p = —2.1 for a power law. This implies that a 
significant fraction of events are caused by MACHOs in the substellar mass range. 
The determined best fitting mass functions are compared with the observed stellar 
mass function as determined by Gould, Bahcall, & Flynn (1996) using the Hubble 
Space Telescope (HST). Both best fitting mass spectra differ by more than 3a levels 
from that of stars in the solar neighborhood. 



2 Models 

We adopt a bar-structured Galactic bulge with a double exponential disk. The Galac- 
tic bulge is modeled by a "revised COBE" model which is based on the COBE model 
(Dwek et al. 1995) except for the central part of the bulge. In the inner ~ 600 pc 
of the bulge, we adopt the high central density Kent (1992) model since the the 
COBE model does not match very well with observations in this region. The disk 
is assumed to have an exponential distribution with vertical and radial scale heights 
of Hz = 325 pc and Hr = 3.5 kpc, respectively (Bahcall 1986). The models are 
discussed in detail by Han & Gould (1995a). 



The transverse velocity is defined as 
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where Vi, Vg, and Vq are the projected velocities of the lens, the source, and the 
observer. In our model, disk MACHOs are assumed to have a Gaussian velocity 
distribution represented by 



f{vi) (X exp -{vi - ^;^)V2(J^ 
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where the coordinates {x, y, z) have their center at the Galactic center and x and z 
axes point to the Sun and the north Galactic pole, respectively. We adopt the disk 



component velocity distribution as -u^^disk = km s~^, cr z^disk = 20 km s~^, -Uy^disk = 
220 km s~^, and Uj^^disk = 30 km s^^. The velocity of the barred bulge model is de- 
duced from the tensor virial theorem (Binney & Tremaine 1987; Han & Gould 1995b) 
with resulting velocity dispersion components {ax', (Jy/, az>) = (115.7, 90.0, 78.6) km s"'' 
Here, the coordinates {x', y' , z') are ahgned along the axes of the triaxial Galactic 
bulge; the longest axis to the x' direction and shortest to the z' direction. The 
projected velocity dispersions are then computed by cr^ = o"^, cos^ 6 + o"^, sin^ 6, 
(Ty — a'^isiv?0 + cr^, cos^ 6', and = (Jz', where 9 ~ 20° is the angle at which 
one views the major axis of the bulge. The resulting values are (cr^;, Uy, Uz) = 
(113.0, 93.0, 78.6) km s~^. We assume a nonrotating bulge, but see Blum (1995). 



3 Mass Spectrum 

We have so far described the models of v and D = DoiDis/Dos, which are the two 
parameters required for the complete determination of the time scale te for a given 
mass m. The next step is to model the mass spectrum of lenses. We test three mass 
functions. In the first model, we assume that log m is Gaussian distributed: 
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where m is expressed in units of Mq. The two free parameters in this mass function 
are the mean and standard deviation of the logarithmic mass: (logm) and a\ogm- 
Note that o"iogm is the width of the mass function not the error estimate in the mass 
range of the mean. We also test the standard power-law mass function, 

/(m) =irm*'e(m-mcut), (3.2) 

where is the Heavyside step function, and K = {p + l)m^~^. Two different 
expressions for the mass function /(m) and /(logm) are related by f{m)dm = 
/(log m)o? logm. In addition, the stellar mass function based on observations with 
the HST (hereafter referred to as the Hubble model) is tested and compared with 
the power-law and Gaussian mass functions. The details of the Hubble model are 
described in §4. 

First, we compute the distribution of time scales f{t'^) for a fixed value of mass 
(e.g., 1 Mq). To find f{t'^), one should weight events by the transverse speed v 
and the cross-section (i.e., Einstein ring radius re). This is because events with 
faster transverse speeds and larger cross-sections are more hkely to occur. Then the 



distribution function of t'^ is computed by 
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where D = Do\D\s/ Dos, and a normalization factor to be discussed below. We 
assume that the bulge is cut off at 4 kpc from the Galactic center (at Rq = S kpc), 
and therefore (imax = 12 kpc. The increase of volume element, and thus increase in 
total number of stars in the volume is assumed to be compensated by the decrease of 
observable stars due to decreasing detectability with distance. This is equivalent to 
the parameter (3 — —1 model of Kiraga &: Paczynski (1994). We take disk self-lensing 
events into consideration by setting the lower limit of source stars to be in equation 
(3.3). 

At the second stage, we compute the actual time scale distribution f{tc) of MA- 
CHOs whose masses are distributed by a function (f){m). The distribution /(te) is 
obtained by convolving two functions, f{t'J and 4>{m): 

fiQ^riiQ dmm}'''ct>{m) dt'J(Q5(t, - t',V^), (3.4) 
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where r]{te) is the detection efficiency, and the factor m^/^ is included to weight events 
by their cross section which is Te oc m^/^. In the computation of /(tc), we assume an 
upper mass limit rriup = 10 Mq. Using the assumed distributions of v, D, and m, 
we find the best fitting parameters of the mass functions; (logm) and aiogm for the 
Gaussian and Wcut and p for the power-law mass spectrum. The best fitting mass 
function is obtained by comparing model and observed time scale distributions. For 
this, we use the maximum likelihood test in which the statistic is computed by 

A^tot 

In L = ^ ln/(te,obs,i), (3.5) 

i=l 

where te,obs is the observed time scale. Prior to the computation of In L, the test distri- 
bution f{te) constructed from equations (3.3) and (3.4) is corrected by the detection 
efficiency r]{te) using the functions provided by Alcock et al. (1995b) for MACHO 
and Zhao, Spergel, & Rich (1994) for OGLE. Since the efficiencies for each group are 
different, the likelihood has been computed separately and summed: 

In L = ^ ^ ln/j(te,obs,i), (3.6) 
j 1=1 



where the subscripts j — 1,2 represent values of MACHO and OGLE, respectively. 
In this test, one determines the best fitting distribution that maximizes In L. The 
likelihood statistic is related to the uncertainty by A(\n L) = a'^/2. Since our con- 
cern is finding best fitting mass spectrum not the total amount of optical depth or 
frequency, we leave the overall normalization ^ as a variable so that the expected 
number of events matches with that of the observed events. 

The distribution of observed time scale is obtained from data measured by MA- 
CHO (Alcock et al. 1995b) and OGLE (Udalski et al. 1994) groups. Out of 55 can- 
didate events detected (45 by MACHO and 12 by OGLE, with 2 overlapping ones), 
we exclude four suspected by MACHO (Alcock et al. 1995b) as being either variables 
or influenced by neighboring variables. For the photometry toward a dense field such 
as Baade's Window, it is expected that blending of stars causes errors in determining 
time scales. The blending is more serious for faint stars and resultant time scales 
determined are shorter than the actual values. To check how seriously the blending 
affects the time scale measurement, we plot the time scale versus V mag of source 
stars in Figure 1. In the figure, data from MACHO and OGLE groups are marked by 
filled and empty circles, respectively. Except three long events detected by MACHO 
group, there does not exist any trend in the figure. We therefore ignore blending in 
our analysis. 



4 Result and Discussion 

The results of the In L computation in parameter space are shown as contour maps 
in Figure 2 for power-law (upper panel) and Gaussian (lower panel) mass functions, 
respectively. In the figure, the best fitting positions in parameter space are marked 
with 'x' and the contours are drawn at la, 2a, and 3a levels; the power-law fit is 
better than the Gaussian by ~ 2.3 a and better than the Hubble by ~ 5.3 a. For 
the Gaussian mass function, we find the best fitting parameters of (logm) = —1.12 
and aiogra = 0.57. The best fitting parameters for the power- law distribution are 
n^'cut = 0.04 Mq and p — —2.1. Therefore, a significant fraction of events might be 
caused by lenses with mass less than the hydrogen-burning mass limit (~ 0.1 Mq) 
regardless of the assumed form of the mass spectrum. 

The determined best fitting mass functions are compared with that of local disk 
stars to check similarities and differences between the two populations. Gould et 
al. (1996) used HST to measure the luminosity function of local Galactic disk stars. 
They then applied mass-luminosity relation from Henry & McCarthy (1990) to obtain 
a local mass function. The Hubble mass function is shown and compared with the 



best fitting power-law and Gaussian mass functions in Figure 3. We liave added the 
wliite dwarf (WD) population to the Hubble model, the bumpy peak at logm ~ —0.3, 
which is modeled based on the observation of McMahan (1989). The mass function 
of WD is normahzed using the mass density pwD ~ 0.005 MqPC~^ (Bahcall 1984). 
Luminous massive stars would have a low scale height and thus only nearby objects 
could contribute to microlcnsing events seen toward Baadc's Window, b ~ —4°. Then 
it is expected one would detect these objects due to both their closeness and high 
luminosity. However, there are no events in which the lens is bright enough to be 
easily detected. Therefore, we cut off the Hubble mass function at the upper hmit of 
~ 2.2 Mq. The best fitting time scale distributions f{te) for the three mass functions 
are shown in Figure 4 and are compared with the observed distribution represented 
by histograms in each figure. Both power-law and Gaussian mass functions fit signif- 
icantly better than the HST model does, implying that the lenses have very different 
mass distribution from stars in the Galactic disk. 

It is puzzling that any type of assumed mass function cannot match very well for 
the long time-scale events (^e > 70 days). Both observed and best fitting f{te) are 
magnified and shown in the smaller boxes inside each panel so that the differences 
between observation and the best fits of each mass functions can be compared more 
easily in the long event region. The expected number of events with tc > 70 days 
are 0.68 and 1.00, for the power-law and Hubble mass functions, respectively. Thus, 
the Poisson probabihties of observing 4 or more long events are ~ 1% for the power- 
law and ~ 2% for the Hubble mass function. The probabihty is even lower (0.55 
events) for the Gaussian. Even though only 4 events have time scale longer than 70 
days, these long events are important because they are responsible for the significant 
fraction of optical depth measured toward the bulge. 

Long time-scale events can be produced under various conditions. First, the long 
events might be caused by a mass function in which a larger fraction of lenses are in 
the higher mass range. If the mass function is a power-law, this can be achieved by a 
higher mcut and lower p. In Figure 5, the expected time-scale distributions for various 
lower values of power (upper panel) and higher mass cut-off (lower panel) arc shown 
and compared to that of best fitting distribution. However, as one or both of the 
conditions are satisfied to explain the long events, the fit deviates seriously in the low 
mass region where the majority of events are located. Therefore, the observed bimodal 
pattern of /obs(^e) is not likely to be due to a lower power or higher mass cut-off than 
our determination. Another explanation would be a MAGHO population with low 
transverse speed. For example, if the velocity dispersion of bulge population is very 
low, a MACHO will cross slowly over the Einstein ring, thus producing very long 
events. However, this explanation does not seem to be plausible because the velocity 
dispersion of the bulge is observationally reasonably constrained, and thus it cannot 



be arbitrarily low. Another possibility is that longer events might be caused if a lens 
is composed of a binary system which is so closely spaced that one would detect it as a 
single lens instead of a binary lens. If the fraction of close binary systems is significant, 
the observed time scale distribution would be significantly different from that expected 
with an assumed single star mass distribution. However, the fraction of binary stars 
whose separations are less than 2 AU, which is a typical Einstein ring radius, is < 25% 
for G type stars in the solar neighborhood (Duquennoy & Mayor 1991). Therefore, in 
the computation we do not take into consideration the modification of the time scale 
distribution due to binary systems. Finally, another explanation would be a bimodal 
mass distribution with the second population composed of heavy dark MACHOs 
(e.g., black holes, neutron stars, or white dwarfs) with low transverse velocity. If, 
for example, there was a kinematically cold population of these objects, they would 
have a low scale height and hence would be observed mostly near the sun for sources 
near Baade's Window, i.e., a few degrees from the Galactic plane. Such a population 
would then have a very low transverse speed. 

It is curious that the best fitting power-law does not deviate much from that of the 
Salpeter mass function, which is valid in the mass range M > 1 Mq. Determination 
of the mass function in the solar neighborhood by Miller & Scalo (1979) indicated 
that the increase in the number of stars becomes shallower (decreasing p) at lower 
mass than the classical Salpeter estimate of p ~ —2.35 for stars with m > 1 Mq. 
Their estimate is p ~ —1.4 in the mass range 0.1 Mq < m < 1 Mq. By contrast, 
our determined power p — —2.1 happen to be close to the classical value, and would 
seem to imply that the classical Salpeter mass function extends to lower mass objects. 
Note, however, that the Hubble function, in which the number of stars decreases as 
mass decreases after passing the maximum at m ~ 0.5 Mq, is inconsistent with such 
a power law. 



5 Conclusion 



We determine the mass spectrum of MACHO detected toward the Galactic bulge 
by modeling velocity and matter distributions. For the best fitting power-law and 
Gaussian mass functions, a significant fraction of events is caused by MAGHOs in the 
substellar mass range. In addition, the mass spectrum seems to differ significantly 
from that of local Galactic disk stars. However, caution is required in interpreting 
these results since the mass spectrum of MACHOs determined here is based on models 
of density and velocity distributions. The parameters of these models, e.g., radial 
and vertical scale height of the disk, total mass of the Galactic bulge, and velocity 
field in the central Galactic bulge, are not unambiguously determined. Different 



models can lead to different conclusion. For illustration, wc arbitrarily set the bulge 
velocity dispersions at high levels of o"^:' = (^y' = cr^' = 100 km s^^. The best fitting 
Gaussian mass spectrum is obtained the same way as described in §3 and the best 
fitting parameters are marked by in Figure 2. However, the differences are small: 
/A(logm) = —0.138 and /AtJiogm = 0.056. An additional uncertainty comes from 
the low detection efficiency for very short events. It is possible that a class of very 
low mass MACHOs is generating very short {t^ < 3 days) events that go completely 
undetected. We have issued this possibility in this paper. The MACHO group is 
presently conducting a "spike analysis" to determine whether such ultra-short events 
in fact occur. 
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Figure 1: Einstein crossing times of lensing events obtained by MACHO and OGLE 
versus V mag of source stars. 





Figure 2: Likelihood contour for the power-law and Gaussian mass functions. The 
contours are drawn at la, 2a, and 3a levels from the best fitting position marked by 
'x'. Th best fitting point with higher bulge velocity dispersion, a — 100 km s""-*^, is 
marked by (see § 5). 




Figure 3: The best fitting mass spectrum for tfie power-law and Gaussian mass 
functions are compared witli tlie Hubble mass function. 
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Figure 4: The respective best fitting time scale distribution functions f{tc) for 
power-law, Gaussian, and Hubble mass functions. They are compared with the ob- 
served distribution, /obs(^e), represented by a histogram. The distribution in the long 
time scale region {te > 70 days) are rescaled and drawn in the smaller boxes in each 
panel. 




Figure 5: The time scale distribution functions for lower value of power p and the 
higher mass cutoff mcut are compared with the best fitting distribution withp = — 2.1 
and mcut = 0.04 Mq. 



